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Abstract: Thermal Diffusion Flowmetry (TDF) (also called Heat Clearance 
Method or Thermal Clearance Method) is a longstanding technique for 
measuring blood flow or blood perfusion in living tissues. Typically, 
temperature transients and/or gradients are induced in a volume of interest 
and the temporal and/or spatial temperature variations which follow are 
measured and used for calculation of the flow. In this work a new method 
for implementing TDF is studied theoretically and experimentally. The heat 
deposition which is required for TDF is implemented photothermally (PT) 
and the measurement of the induced temperature variations is done by 
photoacoustic (PA) thermometry. Both excitation light beams (the PT and 
the PA) are produced by directly modulated 830 nm laser diodes and are 
conveniently delivered to the volume under test by the same optical fiber. 
The method was tested experimentally using a blood-filled phantom vessel 
and the results were compared with a theoretical prediction based on the 
heat and the photoacoustic equations. The fitting of a simplified lumped 
thermal model to the experimental data yielded estimated values of the 
blood velocity at different flow rates. By combining additional optical 
sources at different wavelengths it will be possible to utilize the method for 
non-invasive simultaneous measurement of blood flow and oxygen 
saturation using a single fiber probe. 
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1. Introduction 

Thermal Diffusion Flowmetry (TDF) (also called Heat Clearance Method or Thermal 
Clearance Method) is a longstanding technique for measuring blood flow or blood perfusion 
in living tissues [1-4]. The principles of the methods were first conceived by Gibbs as early as 
1933 [1] and many variations and modifications have been proposed since. The common 
principle of all the variants is that the physical variables blood flow and thermal diffusion in a 
tissue are closely related and the former can be inferred by measurement of the latter. 
Typically, temperature transients and/or gradients are induced in a volume of interest and the 
temporal and/or spatial temperature variations which follow are recorded and analyzed. The 
analysis is usually based on a heat diffusion model which describes the conduction and 
convection of heat in the measured volume. 

Photoacoustic (PA) imaging and spectroscopy is based on measurement of the acoustical 
waves which are generated due to the absorption of modulated light in a tested medium. PA is 
widely known for its excellent performance in vascular imaging [5-7]. The main reason for 
that is the relatively high absorption of hemoglobin in the visible and near-IR spectral regions. 
Photoacoustics enables high resolution imaging at depths where most purely optical imaging 
techniques are not functional due to excessive optical scattering [5]. One very useful attribute 
of PA imaging is its ability to determine the oxygen saturation level of the blood (SO2) [5,6]. 
This is typically implemented by using multispectral PA excitation. 

Recently it was demonstrated that the efficiency of the PA effect strongly depends on 
temperature [8,9] and is consequently affected by variations in flow [10]. In this paper we 
further investigate the PA dependence on temperature and flow and describe a photoacoustic 
version of TDF which is based on this dependence. In contrast with conventional TDF 
methods, in which the temperature is measured via a direct thermal contact with the tissue, we 
utilize the temperature dependence of the efficiency of PA generation for remote monitoring 
of the blood's temperature. Moreover, the same optics which delivers the PA optical excitation 
can be used for photothermal (PT) induction of the heat which is needed for implementing the 
thermal diffusion measurement. This approach offers a number of potential advantages over 
the traditional methods. One advantage is that with no thermal contact and since the light is 
absorbed primarily by hemoglobin, the technique selectively probes the blood and not the 
surrounding tissues. A second advantage is the optional high spatial resolution of the method, 
facilitated by optical as well as acoustical focusing, which enables measurement of the flow in 
a specific blood vessel whereas conventional TDF is typically limited to measurement of the 
global flow in a volume containing multiple vessels. Another very important advantage is the 
ability to use multispectral PA excitation or PT excitation for determining s0 2 simultaneously 
with the blood flow. In this way it will be possible to obtain the metabolic rate of oxygen 
which is an important physical parameter in many clinical situations [11-13]. 

While there were numerous studies of vascular PA imaging, only a few dealt with 
characterization of the blood flow. Most of them were based on the PA Doppler (PAD) effect 
[14-18]. PAD was demonstrated in-vivo in optical resolution mode, when the beam was 
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focused to a size close to a single red blood cell, limiting the penetration depth to less than 
1 mm [14,15]. Implementation of vascular PAD in the acoustical resolution mode, however, is 
rather challenging since the Doppler effect requires heterogeneity of the moving medium. 
When the spatial separation between the red blood cells is small relative to the acoustical 
wavelength, the blood seems homogeneous and the effect is suppressed [14,19]. The method 
described in this paper is not based on the Doppler effect but rather on heat transfer. 
Therefore, it does not require high spatial resolution and can be implemented at relatively 
deep tissues. 

2. Theory 

Thermal Diffusion Flowmetry (TDF) is based on the principle that the tissue temperature 
variations in a Volume Under Test (VUT) during some measurement time, T = r(r,t) 

{ 0 < t < f meas ; r e VUT } , in response to thermal excitation, are closely related to the flow in 

the volume and can be used for its evaluation. There are many ways to induce the thermal 
excitation and to monitor the resulting temperature variations. In this study we propose the use 
of the photothermal (PT) effect for thermal excitation and the photoacoustic (PA) effect for 
monitoring the temperature. The quantitative description of the method is divided into two 
parts. In the first part we will introduce the equations which govern the spatial and temporal 
behavior of r(r,f) in the presence of PT excitation, thermal diffusion and flow-dependent 
heat convection. In the second part we will describe the effect of the photothermally generated 
temperature variations on the PA properties of the medium and how this can be utilized for 
PA monitoring of r(r, /) and consequently of the flow. 

2.1. Photothermal excitation of the VUT 

The temperature in a volume of tissue subjected to photothermal excitation is described by a 
heat equation [20]: 

= V [a(r)Vr(r, t)] - u (r ) V r(r, t) + ^ (l% 0 (1) 
dt - 1 . . J ' > , p(r)C,(r) 

mn/iiirtinn rttnvprtinn ' 



where t is the deviation from the baseline temperature, a is the thermal diffusivity, p is the 
density, C p is the heat capacity at constant pressure, u is the velocity of the blood in the 
tissue and W PT is the PT heat energy deposited in the medium per unit time per unit volume. 
For simplicity, metabolic heat generation was neglected. Generally Jf PT (r, t ) can be written 
as a product of a space-dependent function and a time-dependent function as 
W PT (r,f) = // r PT (r)// ( PT (f) . It is also convenient to use normalization where the area under 

H, PT (f)is 1. 

The temporal variations of r are determined by the three terms on the right hand side 
(RHS) of the equation. The first term describes thermal conduction. The second term 
describes the heat which is added or removed from a given location due to heat convection. 
This term exists only where u * 0 and is closely related to the blood flow. The heat induced 
in the medium due to the absorbed light is described by the third term. 

For an excitation pulse which satisfies the thermal confinement condition [6] (but not 
necessarily the stress confinement condition) Jf PT (r, f) can be expressed as 

W PT (r,f) = Hj T (r)S(t) . In this case the problem can be formulated as an initial value 
problem: 



dT( -*> 1 ) = V [a(r) V r(r,t)] - u (r) V r(r, t) (2) 
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with initial conditions 



r(r,0)= (3) 
p(r)C p (r) 

In general, the solution of Eq. (2) for the initial conditions in Eq. (3) and given boundary 
conditions can be found numerically. Here we concentrate on a simplified lumped solution 
which facilitates the estimation of u from measurement of the temperature. 

2.2. Lumped model 

To gain further physical insight regarding the temporal behavior of r in response to impulse 
excitation we consider a simplified model, in which the VUT is considered a "lump" and its 
thermal state is described by a single parameter [21]. We define the lumped temperature of the 
system as a weighted average over the temperature in the volume: 



f(0 = JJJ/(r)r(r,OrfV (4) 



where /(r) is currently an arbitrary weight function of r, normalized such that its integral 

over the VUT is 1. It is introduced here for use in subsequent paragraphs. Using Eq. (2) we 
obtain 

' A JJJ/(r)V[a(r)Vr(r,0]rfV-JJJ/(r)«(r)Vr(r,r)rfv[f(0 (5) 



dt 

where ¥(r,t) = r(r,t)/f(t) . 

This representation is particularly helpful in situations where the coefficient of f (f ) in the 
RHS of Eq. (5) is independent of time. In this case the decay of f (t) following an impulse 

will be exponential with a time constant ( Conduction + Convection ) . where 

-1 



JJJ/(r)V[a(r)Vr(r,f)]«/V 

(6) 



1 



— - JJJ/(r)«(r)Vf(r,0^ 

VUT 

An important observation that will be used to analyze the experimental results is that in the 
case of uniform flow f COIlv( , ction becomes inversely proportional to the velocity. In many 
practical situations a lumped model is used even if the characteristic decay times are not 
strictly constants. The suitability of the model depends on the thermal and geometrical 
properties of the VUT and its accuracy can be evaluated empirically. The usefulness of the 
lumped model in formalizing the thermal dynamics of a blood vessel phantom and its 
potential role in flow measurement is demonstrated in Section 4. 

A simple but compelling description of the lumped model can be made in terms of an RC 
circuit (Fig. 1). In these terms the illuminated volume is conceived as a thermal "capacitor" 
C themiai which is charged by the PT source. The absorbed energy is dissipated through two 

parallel resistors: /? conducUon and R comection ■ The temporal behavior of the thermal circuit is 
determined by the time constant f cff = ( faction +Cv=ction)" 1 . where f conduction = C thcrmal fi condllctlon 
and Cvec™ = Qcrmai Convection • In the case of zero flow (stationary blood) the convection 
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resistance is infinite and all the energy is dissipated through /? conduction ■ As the flow increases, 
the convection resistance decreases and more energy is dissipated through /? convcction , until, for 
very high flows, it completely shortens the capacitor and prohibits its charging. Moreover, in 
light of the definition of f colldllction in Eq. (6) it can be seen that if the velocity is uniform the 
convection resistance is inversely proportional to the velocity and can be expressed as 
R R \u\ 

convection convection / | | ' 

reconduction Rconvcction 

\ / 




Cthcrmal 

Fig. 1. RC circuit description of the lumped model. 

2.3. Photoacoustic monitoring of the temperature variations 

The generation and propagation of pressure waves in a fluid due to the photoacoustic effect, in 
the so-called thermal-confinement regime, can be described by the following inhomogeneous 
wave equation [22]: 

^4^-cVVr,0=r^M (7) 
dt 2 dt 

where p(rj) is the pressure field, c is the velocity of sound, T = fie 1 C p is the Griineisen 
coefficient, p is the volumetric thermal expansion coefficient, C p is the specific heat and 
Jf FA (r,t) is the PA heat energy deposited in the medium per unit time per unit volume. Here 
as well the heat rate, 9f FA (r,t) , can be written as a product of a space-dependent function and 
a time-dependent function as Jf PA (r,t) = H^ A (r)H^ A (t) . Again, the area under //, PA (f) is 
normalized to 1 . 

To use PA for temperature monitoring we take advantage of the temperature dependence 
of r which, for limited temperature deviations, can be approximated by [8,9] 

r(T)*r(T 0 ) + b(T-T 0 )=T 0 +bz (8) 

where T is the temperature, T 0 is the baseline temperature, r is the deviation from T a , r o is 
the Griineisen coefficient at T 0 and b is the slope of the dependence curve. As described in 
Section 4, we experimentally found that in blood the relative change in F, 
[Ar(r)/r„] It = b/r 0 , is approximately 3.5% per degree Celsius in the range 

25° < T < 35° . This is similar to the value measured in water and is attributed mainly to the 
increase in the thermal expansion coefficient [20,23]. 

In the present form Eq. (7) describes the situation where T is constant. In the presence of 
PT excitation F becomes a function of both space and time due to the spatial and temporal 
dependence of r(r, t) . Since T determines the fraction of deposited heat energy which is 
converted to pressure at each point in space at a given time, the Poisson solution to Eq. (7) in 
the case of an optical impulse at t = t ' becomes 
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p(r,t-t',t') 



1 d 
4xc dt 



(9) 



A(t-f) 



Here A(t) denotes a surface where |r-r'| = ct . Note that the Poisson solution at a given 

position, p(r 0 ,t — t',t') , is no longer a function of t — t' alone but rather it depends on the time 

of excitation t' and 7/ r PA (r) is weighted according to the spatial dependence of T. It is also 

assumed that T changes very slowly with respect to the characteristic PA response time so it 
can be assumed constant during the response. For an excitation waveform other than an 
impulse the pressure field can be obtained by 



p(r,t)= J Hj A (t')p(r,t-t',t')dt' 



(10) 



2.4. Temperature monitoring using sinusoidal PA excitation 
Consider the case where the PA excitation is sinusoidal, namely, 



H ,'A W = J_[i + C0B(4DtAf )] = 1 [ 2 



+ e r 



+ e 



(11) 



where a> PA is the radial frequency of the modulation and f total is a normalization constant 
which takes into account the fact that in reality the excitation has a finite duration. 
Substituting in Eq. (10) and retaining only the term centered at +« PA yields 



2t 



p(r,f,t-f)dt" 



(12) 



To proceed, we substitute the exact form of p(r,t",t — t") from Eq. (9) and the expression 

for T from Eq. (8) and arrive after several mathematical steps to the following simple result 
(see Appendix A for details): 



i- 



-m 



(13) 



where we assumed that r o is spatially uniform, the observation point is at the origin, Ae' 
is the response obtained from Eq. (12) in the absence of photothermal excitation and 

r(t) = jjjf(.r)T(t,r)dV (14) 



where 



H PA ( r ) <,-■*■** W« 



/(r) 



jjj7(r)JV 



(15) 



It can be seen that the PT excitation leads to Amplitude Modulation (AM) of the otherwise 
purely sinusoidal PA response. The carrier amplitude, A, can be found by turning off the PT 
excitation. With A and the relative change in the Griineisen coefficient known, the 
temperature modulation term can be found from the measured PA response according to 



r 

V 1 o J 



— 1 



(16) 
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For an impulse PT excitation, in the cases where the lumped model applies, fit) takes an 
exponential form according to Eqs. (5) and (6): 



fit) = f (0) exp 



f P 



V f efl J 



(17) 



The Fourier transform of Eq. (17) yields the Modulation Frequency Response (MFR): 

f(a>)=f(0) hs (18) 

It can be seen that in the lumped model approximation the photothermal MFR takes the 
form of a complex Lorentzian. The applicability of this result to our experimental conditions 
is demonstrated in Section 4. By fitting the measured MFR to a Lorentzian, f eff was found 
and used for estimation of the velocity. 

3. Experimental setup and measurements 

The experimental setup (Fig. 2) comprised a pair of 830 nm Laser Diodes (LD's), each with a 
maximum peak power of 200 mW. The outputs of the LD's were polarization-combined and 
the resulting beam was coupled to a multimode fiber. The intensities of both LD's were 
sinusoidally modulated by direct modulation. One LD was modulated at a frequency of 
f PA = 0.5 MHz to generate a PA signal. To characterize the MFR the other LD was 

sinusoidally modulated at discrete frequencies in the range of f PT e [0.1 ,20] Hz . This choice 

of a frequency-domain characterization of the modulation response rather than a time-domain 
approach based on a thermal impulse was made to enable a higher signal to noise ratio. The 
unconnected end of the fiber was brought in close proximity to a vessel phantom comprising a 
transparent Tygon tube (1mm inner diameter) filled with defibrinated sheep blood. The size of 
the illuminated volume was roughly 0.75 mm 3 . The central part of the tube was immersed in 
water while one end was attached to a syringe pump and the other end to a container. 

Detection of the PA signal was performed using a non-focused ultrasonic immersion 
transducer with a center frequency of 0.5 MHz (Olympus IR-0008-P) followed by a low-noise 
preamplifier (Olympus 5660C). The amplified signal was sampled using a digital acquisition 
card and recorded by a PC. The sampling rate was 1.25 MS/s and the acquisition time was 20s 
per data set leading to spectral resolution of 0.05 Hz. Each measurement of PA response was 
repeated twice and the average response was calculated. A Fourier transform of the PA 
response yielded a peak at f PA and modulation sidebands at f PA ± f PT . 

As a preliminary stage, the temperature dependence of the PA response in stationary blood 
was characterized. This was performed by a simultaneous measurement of the blood's 
temperature and its PA response. To measure the temperature, a thermocouple was inserted 
into the tube and brought in close proximity to the illuminated volume. The position of the 
thermocouple was determined as a compromise between the need to place it as close as 
possible to the VUT center and the attempt to avoid PA generation by its direct illumination. 
To minimize the inaccuracy, the readings were taken only after sufficient time when the 
transient component in the temperature signal decayed and the system reached a steady state. 
It should be noted that implementation of the PA-TDF method does not require a 
thermocouple nor does it use the value of the Griineisen temperature dependence that was 
estimated from the thermocouple measurements. 

To experimentally test the PA-TDF method, the PA modulation frequency response was 
measured for several flow rates by changing the PT modulation frequency in the range of 

f PT e [0. 1 , 20] Hz . The amplitudes of the two PT sidebands were averaged for each value of 

/ PT to generate the modulation frequency response (MFR). This was repeated for both 

stationary blood and flowing blood at 5 different flow rates between 0.05 and 1 ml/min. The 
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corresponding blood velocities were between 1.06 and 21.2 mm/sec for the tube that was 
used. 

For each blood velocity u, a normalized Lorentzian was fitted to the measured MFR 
according to Eq. (18). This procedure yielded the effective time constant as a function of 

velocity, f cff («) . The conduction characteristic time, f conductio „, was obtained from the 



stationary measurement f conductlon 
calculated according to t c 



■■ f cff (u = 0) . With f eff and f ( 



convection 



co„duc,o„ known, f convcction could be 
= (tea ~ Conduction ) • As described in the next section, in our 
tested phantom ? convection was found to be inversely proportional to u . The parameter of 

proportion, having units of length, can be interpreted as a measure to the spatial dimension of 
the VUT in the direction of the flow. 



830 nm Directly 
Modulated Laser 
Diodes (x 2) 



Arbitral 
Wavefoi 
Generad 



rm 
tor 




Fig. 2. The experimental setup. 



4. Experimental results 

The direct measurement of temperature using the intra-tube thermocouple enabled 
characterization of the temperature-dependence of the PA response. Temperature variations 
were induced in the measurement volume by sinusoidal PT modulation at a frequency of 
0.1 Hz. To scan a wide temperature range the peak power of the photothermal-LD was tuned 
between 0 and 200 mW in increments of 25 mW. At each level of peak power the readings of 
the thermocouple, after a transient time, oscillated in correspondence with the PT optical 
excitation. We denote the amplitude of these temperature oscillations as r AC [C°] and the 
deviation of their mean value from the baseline temperature (the thermocouple output in the 
absence of photothermal excitation) as r DC [C°]. In parallel with the PT excitation, the 

volume was illuminated by the photoacoustic-LD whose power was sinusoidally modulated at 
a frequency of 0.5 MHz. The presence of PT excitation led to modulation of the generated PA 
signal. The recorded PA signal was spectrally analyzed and its spectral peaks were identified 
and recorded. An example of the PA response spectrum in the presence of PT excitation at 
frequency of 2 Hz is plotted in Fig. 3. The modulation sidebands due to the PT excitation at 
/ PA ± are evident. 

The dependence of the central PA peak on r DC normalized by its value without thermal 
excitation, is plotted in Fig. 4(a). An estimate to this ratio can be found from Eq. (13) by 
approximating f(?)«r DC (where x(f) denotes the time average of x(?)). It can then be 
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-8 -6 -4 -2 0 2 4 6 8 10 
Frequency detuning [Hz] 

Fig. 3. The spectrum of a PA response with PT modulation at 2Hz. 

shown that the graph in Fig. 4(a) should follow 1 + (b/T 0 ) r DC . The linear fit in Fig. 4(a) 

implies a slope of b/T 0 = 0.034 [1/°C] . This result, which manifests a temperature 
dependence of 3.4% per degree Celsius of the Griineisen coefficient of blood, is in close 
agreement with the value obtained for water which was around 3.85% [1/C°] [20,23]. The 
dependence of the modulation peak on r AC together with a linear fit is plotted in Fig. 4(b). In 

this case the anticipated slope is (&/r o )/2 since the value is taken from a modulation 

sideband whose amplitude is half the amplitude of the modulating sinusoidal signal. The fitted 
slope, however, was 0.015 which is slightly lower than the expected value. The deviation is 
attributed to the very long response time of the thermocouple, which leads to underestimation 
of t ac even in the relatively low modulation frequency of 0.1 Hz. The thermal properties of 
the VUT, as well as the excitation frequency, affect the ratio between r AC and r DC , and as 
will be seen in the next paragraphs, this ratio should increase with the increase of flow. 

Normalized modulation frequency responses for stationary blood (q 0 =0 ml/min) and 5 

flow rates of q 15 e [0.05, 0.1, 0.2, 0.5, l] are shown in Fig. 5. For the tube that was used these 
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Fig. 4. Simultaneous measurements of temperature and PA response at different PT modulation 
power: (a) PA amplitude of carrier frequency, normalized by amplitude without PT 
modulation, vs. the average temperature r D c- Experimental data marked by squares, linear fit in 
dashed black, (b) The amplitude of 0.1 Hz modulation peak, normalized by the carrier 
amplitude without PT modulation, vs. rxc- 
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flow rates corresponded to the following blood velocities: w o s e [0, 1.06, 2.1, 4.2, 10.6, 21.2] 
[mm/s]. All responses were normalized so that their value at / PT = 0.1 Hz is 1 and they are 

given in dB scale. In the lumped model approximation the MFR's follows a Lorentzian 
behavior. In accord with Eq. (18) the normalized responses were fitted to 

l + (2^0.1f eff ) 2 J 1/2 y^l + (2;r/p T f eff ) 2 J 1/2 . The suitability of the lumped model for the 

experimental phantom can be evaluated from the Lorentzian fits in Fig. 5. The agreement is 
generally good for the four lowest velocities but rather large deviations can be observed in the 
plots of the two highest velocities. These deviations are attributed to the distribution of the 
fluid velocity in the tube. As the flow rates increases, the velocity profile, which follows a 
parabolic shape in a laminar flow regime, becomes more and more peaked and the 
approximation of uniform velocity, which is implicitly assumed in the lumped model, 
becomes increasingly less accurate. Clearly, to describe the MFR for non-uniform flow the 
lumped model should be modified. It is shown below, however, that despite the discrepancy in 
the MFR, the estimated velocities for non-uniform flow fit well with the expected mean 
velocities. 
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-5 
-10 

I " 15 
o 

■Si 

o -20 
a, 

i/i 
a 

< -25 
-30 
-35 
-40 

0 2 4 6 8 10 12 14 16 18 20 
Modulation Frequency [Hz] 

Fig. 5. Normalized PT-PA modulation frequency responses for stationary (blue) and velocities 
of 1.06 (pink), 2.1 (green), 4.2 (red), 10.6 (black) and 21.2 (gray) mm/s. Solid lines with 
markers are experimental data. Dotted lines are Lorentzian fits. 

The best fit time-constants were, as expected, a decreasing function of the blood velocity: 
# e [0.838, 0.405, 0.208, 0.129, 0.054, 0.029] [s]. Using f convection =(£ -C^j' (where 
'c„» d uc„o„ - feft (« =0)) yielded Cec™ e K 0-784, 0.277, 0.153, 0.058, 0.03] [s]. As can be 
seen in Fig. 6, these results manifested a linear relation between t^ mcctjon and the fluid's 
velocity. The optimal proportion parameter, which we denote as Z VUT , was found by 
minimizing the error between the "real velocity," u, and the "estimated velocity," 
"e St , m a,ed = 'vut Aco„vect,on ■ The error was defined as ^["^vmALve.*™,] 2 (k = 0..5). The 
physical interpretation of l vm (and the origin of its notation) is obtained by recalling that 

'convection ^ s trie time it takes for the temperature to return to its steady state value, after a short 
heating pulse, in the absence of heat conduction. In this case the only heat-clearance route is 
flow and a steady state is reached as soon as the heated fluid leaves the VUT. Hence, 
/ VUT « u ■ f convcction can be interpreted as a length parameter describing the dimension of the 
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VUT in the direction of the flow. The value of Z VUT which minimized the error between 

"estimated u (f° r au> measured velocities) was found to be 0.62 mm,, which is in general 
agreement with the Full Width Half Max (FWHM) of the optical beam. It should be noted, 
however, that the VUT does not have sharp boundaries and analytical determination of / VUT 
from the spatial distribution of the PT excitation requires additional analysis. Interestingly, it 
can be seen that the linear dependence extends even to velocities where the Lorentzian fit 
showed relatively large deviations from the measured frequency responses. The maximum 
error between the estimated and the real velocity was 0.3 mm/sec. Also plotted in Fig. 6 
(green circles) are velocities that were estimated while ignoring thermal conduction, namely, 

using w cstimatcd =/vtrrAeff w ^ tri tne best-fit Z VUT for this case. As can be seen in Fig. 6 this 
resulted in increased errors, in particular in the case of low velocities. The maximum error in 
this case was 0.75 mm/sec. 

Another graph in Fig. 6 (blue crosses) shows the estimated velocity based on data from 
two PT modulation frequencies only (1 and 15 Hz). This can lead to a significant reduction in 
measurement time (<2 s) at the expense of increased estimation error (in particular in the high 
velocities). 




10 15 
Real velocity [mm/s] 

Fig. 6. Estimated velocity ( = /' tollItcll01 , ) vs. real velocity: estimation based on the complete 

lumped model (red circles). Estimation based on a lumped model which ignores thermal 
conduction (green circles). Estimation based on only 2 PT modulation frequencies: 1 Hz and 15 
Hz (blue x). Dotted black line indicates the 45° line. Inset: zoom-in on the range of low 
velocities. 

The system's sensitivity to temperature modulation can be estimated based on the 
measured value of b/T 0 and the photoacoustic SNR at the carrier frequency. Using Eq. (13) 
the ratio between the PA response at the modulation sideband / PA ± / M and the PA response 
at the earner frequency / PA can be approximated as 



P(/pa±/pt) {blT 0 )r AC l2 



P(/pa) i + (^r 0 )r D 



:(fc/r 0 )r AC /2 



(19) 



where it was assumed that (Z?/r o )r DC « 1 . Therefore, the minimum detectable temperature 
oscillation amplitude can be estimated by 
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with SNR(f PA ) defined as the ratio between p(f FA ) to the noise power at / PA ± / PT . In Fig. 3 

for example, SNR(f PA ) « 1MB , which leads to r AC min «0.01°C. Note, however, that this 

sensitivity was obtained when the PA excitation was at maximum power. Any decrease in the 
PA excitation power will result in an increase in the minimum detectable temperature 
oscillations. 

Finally, the system's ability to operate at lower input powers was tested. For PT 
modulation frequency of 0.1 Hz, PT power of 40 mW, PA power of 30 mW and stationary 
blood, a PT modulation peak of ~5dB above noise level was obtained. These levels give a 
rough lower bound to the required optical power. It corresponded to t dc of 0.9°C (out of 

which about 0.5°C was contributed by the PT excitation) and to t ac of about 0.1C°. As 

described in Section 3, the readings of the thermocouple possessed some inherent inaccuracy. 
By waiting sufficient time before taking a reading (until the transient component in the 
temperature signal decayed) and carefully positioning the thermocouple these inaccuracies 
were minimized. With these precautions it is estimated that the total increase in temperature 
did not exceed 2°C. Such temperature increments are comparable with the ones which are 
induced in other TDF methods, and are well below the threshold of thermal damage [24]. 
Moreover, the proposed measurement technique can potentially be performed using pulsed 
PA excitation, provided that the pulse repetition frequency is significantly higher than the 
photothermal modulation frequency. In this case it is expected that a much lower mean PA 
input power will be required. 

5. Discussion 

The experimental results described above provided a first proof of concept to Photoacoustic 
Thermal Diffusion Flowmetry. They demonstrated, albeit in-vitro, the potential of the method 
to become a viable tool in vascular characterization. The main advantage of the proposed 
method over PA Doppler flowmetry is its ability to operate in the acoustic-resolution regime 
(rather than in the optical-resolution photoacoustic-microscopy mode), which is challenging 
for PAD due to the homogeneity of the blood. The ability to operate in the acoustic-resolution 
regime leads to a much bigger measurement depth and simplifies the implementation of the 
measurement. 

It was found that reasonable estimation of the blood velocity, in particular in the case of 
relatively high velocities, can potentially be obtained without a need for calibration. The only 
needed parameter is / VUT . While this parameter was obtained from a minimization process in 
our experiment, it may be possible to estimate it based on the dimensions of the optical beam 
and the estimation of the optical and thermal propagation in the tissue. For low velocities, 
where heat conduction plays a significant role, more accurate results can be obtained if the 
conduction time constant is used to correct the otherwise overestimated velocities. The 
conduction time can be found via a reference measurement of stationary blood. In-vivo, this 
can be implemented by a temporary occlusion of the vessel. 

Another practical consideration is the duration of the measurement. It is of course desired 
to implement the measurement in the shortest time possible. While the experimental 
procedure that was used included measurements in multiple PT excitation frequencies, it was 
found that reasonable estimation of the velocity can be obtained from measurements at two 
frequencies only. There are three main factors that will affect the choice of the PT excitation 
frequencies in realistic scenarios: the measurement duration, the sideband amplitude and 
undesired spectral broadening of the PA peak at the carrier frequency. The latter effect is 
expected to be generated by PA-induced temperature variations, flow variations in the VUT or 
patient movement during the measurement. While in the present experiment the broadening of 
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the carrier was not significant, it is expected to be more pronounced in in-vivo measurements 
and to limit the minimum useable PT frequency. 

Our study was focused on the case of a single vessel. The same approach, however, can be 
applied to a tissue containing a large number of blood vessels with intricate structure. In this 
case the underlying model and equation can be replaced with a bio-heat model [3,4]. This case 
is more similar to conventional TDF techniques, which generally measure a global parameter 
such as the tissue perfusion. Hence, it is to be expected that the proposed measurement 
approach could also be implemented over a volume containing multiple vessels and make use 
of the excellent analytical tools that were developed to make TDF a viable quantitative 
measurement method of blood flow [2—4]. 

Another aspect that must be taken into consideration is the safety of the method. Here as 
well, the similarity to TDF can be used for assessment of the possibility for safe 
implementation. TDF implementations which heat the tissue for flow measurement are 
already commercially available [25]. Accordingly, it is reasonable to assume that the proposed 
method can be implemented under similar restrictions. Experimentally we observed that it was 
possible to obtain good SNR with heating of less than 2C°, which is well below the threshold 
for thermal damage [24], using sinusoidal modulation of both diodes. Although an additional 
source for the PA excitation is required, as mentioned in Section 4, the CW modulation of the 
PA source is not essential to the method and can be replaced by a pulsed source that produces 
considerably less heating. The pulse repetition rate (PRF) of such PA source needs to be high 
enough compared to the thermal clearance rates (i.e. >100 Hz). Spectrally, the PA signal will 
show discrete harmonics separated by the PRF while the PT modulation sidebands will be in 
the near vicinity of each PA harmonic. Further optimization of the system, including focusing 
of the optical beam and/or ultrasonic transducer, may also improve the sensitivity and allow 
detection at lower optical powers. 

6. Conclusions 

In this paper a new method for implementing Thermal Diffusion Flowmetry was proposed, 
analyzed and tested experimentally. The method is based on photothermal excitation of a 
VUT and photoacoustic measurement of the resulting temperature variations. The PA-TDF 
method enabled measurement of blood velocity in a phantom vessel in a velocity range of 0- 
22 mm/s with errors smaller than 0.3 mm/s for all velocities. The relation between the velocity 
and the measured PA signal was derived from the heat and the PA equations via the use of a 
simplified lumped model. As the method is based on light absorption its multi-spectral 
implementation would provide significant spectroscopic information. In this way, for 
example, it is expected to accomplish a simultaneous measurement of both blood velocity and 
its oxygenation level (s0 2 ). 

Appendix A: PA response with sinusoidal PA excitation — detailed derivation 

In this appendix the response to sinusoidal PA excitation in the presence of a general PT 
excitation (Eq. (13)) is derived. As described in Subsection 2.4, the PA excitation is taken as 

Hf A (f ) = — [l+ cos (<y pA f )] = —^[2 + e"^' + ] (A. 1) 

^total ^total 

where ca PA is the radial frequency of the PA modulation and f total is a normalization factor 
which stems from the finite duration of the excitation. Substituting in Eq. (10) and retaining 
only the term centered at +ffl PA the PA response becomes 

■"total -oo 

Using the Poisson solution from Eq. (9) yields 
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The temporal derivative in Eq. (A. 3) can be eliminated via integration by parts: 



4;rc2f„ 



J e W | 



r(r',f-f")//; A (r') 



r-r 



(A.4) 



'total -oo ^('") 

Since T varies very slowly with respect to the typical PA impulse response T(r',t-t") 
can be replaced with T{r',t). Equation (A.4) comprises two nested integrals. The outer 
integral sums over t" . The inner integral sums over a surface of a sphere with a time 
dependent radius |r -r'| = c?" . Hence, each point in time corresponds to a different surface 
and the time integration can be transformed to the volume integral: 



jco PA e ia ™' m r(r',t)H* A (r')e 



-dV 



r-r 



(A.5) 



'total VUT 

Now substituting the form of T given by Eq. (8), and assuming r o is spatially 
independent, 



r-r 



'total VUT 

Setting, without loss of generality, the observation point at the origin (r = 0), the PA 
response becomes 



P^(0 = A 



l +—f(t) 

r 

1 n 



,i<"PA' 



where 



: (t) = $f(r)T(t,r)dV 



f(r> 



f(r) 



(A.7) 

(A.8a) 
(A.8b) 



(A. 8c) 
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